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Abstract 


A  comprehensive  set  of  equations  is  derived  for  correlations  and  spectra 
of  atmospheric  turbulence,  useful  for  wave  propagation  calculations  and 
other  applications.  Three  basic  turbulence  models  are  considered:  the  Gaus¬ 
sian,  von  Karman,  and  Kolmogorov  models.  Two  extended  forms  of  the 
von  Karman  model  are  also  described:  the  Mann  model  for  a  shear-driven 
atmospheric  surface  layer,  and  the  Hunt/ Graham/Wilson  model  for  a  con¬ 
vective  boundary  layer.  A  new  method  for  synthesizing  random  fields  from 
an  inhomogeneous  spectral  model,  called  the  generalized  random-phase 
method,  is  described  and  applied  to  the  various  turbulence  models. 
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1.  Introduction 


Atmospheric  turbulence  interferes  with  the  coherent  propagation  of  acous¬ 
tic  and  electromagnetic  waves.  As  a  result,  turbulence  becomes  important 
in  many  problems  involving  detection,  identification,  and  location  of  tar¬ 
gets.  Statistical  models  of  turbulence,  which  are  useful  for  incorporation 
into  wave  propagation  and  scattering  calculations,  are  the  subject  of  this 
report. 

The  report  is  intended  mainly  as  a  reference  document:  it  brings  together 
a  comprehensive  collection  of  results  from  turbulence  spectral  modeling 
relevant  to  wave  propagation  calculations.  Many  of  the  results  have  not 
been  previously  published,  or  are  scattered  among  diverse  publications. 
A  secondary  purpose  of  the  report  is  to  provide  a  unified  and  systematic 
derivation  of  several  statistical  turbulence  models,  as  well  as  methods  for 
synthesizing  random  fields  from  these  models.  The  report  should  be  useful 
to  readers  having  only  modest  familiarity  with  turbulence  theory  (for  ex¬ 
ample,  readers  whose  specialty  is  wave  propagation)  by  providing  a  con¬ 
densed  tutorial  on  turbulence  spectral  theory  and  modeling. 

The  motivation  for  assembling  these  results  is  the  strong  interest  that  has 
developed  during  the  past  several  years,  particularly  in  acoustics,  in  incor¬ 
porating  realistic  turbulence  models  into  scattering  calculations.  This  inter¬ 
est  has  been  driven  by  the  developing  capability  of  computers  and  software 
to  accurately  calculate  atmospheric  turbulence  effects  on  sound  fields.  An 
important  example  application  is  the  scattering  of  sound  into  an  acoustic 
shadow  zone.*  Shadow  zones  are  created  by  topographic  features  (such  as 
hills)  or  by  refraction.  Recent  research  (Gilbert  et  al.,  1990,  Juve  et  al.,  1994, 
Gilbert  et  al.,  1996)  has  demonstrated  that  most  sound  energy  reaching  a  re¬ 
ceiver  in  an  acoustic  shadow  is  scattered  there  by  turbulence,  and  that  the 
simple,  Gaussian  models  favored  in  the  past  do  not  describe  the  scattered 
signal  levels  well.  Accurate  calculation  of  sound  fields  therefore  requires 
realistic  turbulence  models. 

By  turbulence  model,  I  mean  a  set  of  equations  or  an  algorithm  that  mimics 
statistical  or  dynamical  properties  of  turbulence.  Because  of  the  complex 
nature  of  turbulence,  there  is  no  perfect  model.  Therefore,  the  usual  ap¬ 
proach  to  studying  turbulence  is  to  try  to  develop  approximate  models  for 

*In  analogy  to  an  optical  shadow,  an  acoustic  shadow  zone  is  defined  as  a  region  not  di¬ 
rectly  penetrated  by  rays  of  sound. 
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the  particular  properties  of  interest.  This  report  is  limited  to  second-order, 
two-point  statistics  of  turbulence:  correlation  functions,  spectra,  and  struc¬ 
ture  functions.  These  are  the  statistics  normally  needed  to  calculate  mean 
sound  levels  in  acoustic  shadow  zones,  and  to  calculate  signal  coherence 
at  an  array  of  sensors.  Some  applications,  such  as  the  modeling  of  sound 
level  fluctuations  in  shadow  zones,  require  turbulence  statistics  higher  than 
second  order;  however,  such  higher  order  models  are  beyond  the  scope  of 
this  report. 

In  this  report  I  develop  second-order  statistical  model  equations  by  using 
a  prescribed  equation  for  the  specific  energy  spectrum*  as  a  starting  point. 
I  chose  this  approach  because  it  leads  to  well-determined  and  straightfor¬ 
ward  modeling  procedures.  Alternative  procedures  for  developing  turbu¬ 
lence  models  are  however  possible:  for  example,  one  can  define  models  by 
starting  from  correlation  functions.  A  disadvantage  with  this  alternative  is 
that  it  can  lead  to  confusion  when  vector  fields  are  involved.  Furthermore, 
I  show  in  this  report  how  the  energy  spectral  approach  can  save  some  work 
in  deriving  model  results  needed  for  scattering,  particularly  for  an  impor¬ 
tant  quantity  called  the  two-dimensional  (2D)  correlation  function. 

Regardless  of  how  the  turbulence  model  is  defined,  the  fundamental  diffi¬ 
culty  in  modeling  second-order  turbulence  statistics  lies  in  realistic  model¬ 
ing  of  the  energy  subrange.  This  part  of  the  turbulence  spectrum,  consisting 
of  the  largest  eddies  in  the  flow,  is  directly  affected  by  the  instability  driv¬ 
ing  the  turbulence.  Since  instability  mechanisms  that  drive  turbulence  (pri¬ 
marily  wind  shear  and  density  gradients)  are  inherently  directional  forces, 
statistics  of  the  energy  subrange  eddies  are  anisotropic  (dependent  on  di¬ 
rection).  Statistics  of  the  energy  subrange  are  also  highly  flow  dependent. 
The  anisotropy  and  nonuniversality  hamper  development  of  universal  pa- 
rameterizations  for  energy  subrange  statistics. 

Eventually  the  energy  subrange  eddies  break  down  into  the  smaller  eddies 
of  the  inertial  subrange,  and  the  "memory"  of  the  instability  mechanism  that 
created  the  turbulence  is  lost.  As  a  result,  the  inertial  subrange  is  isotropic 
and  amenable  to  universal  parameterizations.  Using  Kolmogorov's  (1941) 
scaling  hypotheses,  it  is  not  difficult  to  develop  realistic  turbulence  models 
for  the  inertial  subrange. 

The  inertial  subrange  eddies  continue  to  break  down  into  smaller  ones,  un¬ 
til  they  are  so  small  that  molecular  dissipation  processes  convert  the  tur¬ 
bulent  energy  to  heat.  This  transition  occurs  at  eddy  sizes  of  1  to  10  mm 
in  the  atmosphere.  These  small,  dissipating  eddies  make  up  the  dissipation 

*  Specific  here  means  per  unit  mass.  For  example,  the  specific  kinetic  energy  of  a  moving 
parcel  of  fluid  is  u^/2,  where  v  is  fhe  speed  of  the  parcel. 
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subrange.  Modeling  of  the  dissipation  subrange  is  not  considered  in  this  re¬ 
port,  because  sound  waves  that  are  useful  for  long-range  detection  have 
characteristic  wavelengths  that  are  much  larger  than  the  dissipation  sub¬ 
range  eddies  and  therefore  play  an  insignificant  role  in  scattering. 

I  consider  several  turbulence  models,  which  capture  the  turbulence  fea¬ 
tures  discussed  above  with  varying  degrees  of  realism: 

1.  Isotropic  Gaussian  model.  This  is  the  simplest  reasonable  model  from 
an  analytical  standpoint.  All  the  desired  results  are  easily  obtained 
in  closed  form.  The  disadvantage  of  the  Gaussian  model  is  that  it  is 
unrealistic  for  the  inertial  subrange.  The  energy  in  high-Reynolds- 
number  turbulence,  such  as  the  atmosphere,  is  distributed  over  a 
broad  range  of  spatial  scales,  although  in  the  Gaussian  model  all  the 
energy  is  concentrated  at  a  single  length  scale.  The  Gaussian  model  is 
discussed  in  this  report  mainly  for  reference  purposes.* 

2.  Von  Karman  model.  This  model  has  a  realistic  inertial  subrange.  The 
energy  subrange,  although  mathematically  well  behaved,  is  some¬ 
what  unrealistic  since  it  is  isotropic.  All  second-order  statistics  of  in¬ 
terest  can  be  obtained  analytically,  although  the  mathematical  manip¬ 
ulations  become  rather  complicated. 

3.  Kolmogorov  model.  This  model  is  intended  only  for  the  inertial  sub¬ 
range,  for  which  it  is  both  realistic  and  simple.  One  drawback,  in 
comparison  to  the  von  Karman  model,  is  that  many  of  the  equa¬ 
tions  diverge  in  the  energy  subrange.  Therefore  some  care  must  be 
taken  in  applying  the  Kolmogorov  model.  All  results  for  the  Kol¬ 
mogorov  model  actually  can  be  obtained  from  calculations  of  high- 
wavenumber  approximations  to  the  von  Karman  model,  as  well  as 
by  other  methods.  The  Kolmogorov  model  can  be  expressed  entirely 
in  terms  of  a  single  parameter,  called  the  structure-function  parame¬ 
ter.  This  feature  has  proven  convenient  in  previous  experimental  and 
theoretical  studies  of  scattering. 

4.  HGW  (Hunt/Graham/Wilson)  model.  Hunt  and  Graham  (1978, 
1984)  developed  a  method  to  determine  the  modification  of  free- 
space  spectra  due  to  blocking  (damping  of  the  vertical  velocity)  at 
a  plane  boundary  such  as  the  ground.  The  boundary  was  modeled 
as  a  flat,  free-slip  surface.  In  previous  work  (Wilson,  1997c),  I  applied 
the  Hunt  and  Graham  method  (with  the  von  Karman  model  used  for 

*One  advantage  of  the  Gaussian  model  is  that  it  is  particularly  easy  to  formulate  an 
anisotropic  version.  As  a  result,  it  may  work  reasonably  well  for  the  energy  subrange  in 
some  situations  (Wilson  &  Thomson,  1994).  In  this  report,  however,  I  consider  only  the 
isotropic  Gaussian  model. 
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the  free-space  spectra)  to  the  atmospheric,  convective  boundary  layer 
(CBL).  The  HGW  model  is  an  example  of  a  vertically  inhomogeneous 
(height-dependent)  model.  A  disadvantage  of  the  HGW  model  is  that 
many  statistics  of  interest  must  be  obtained  by  numerical  integration. 

5.  Mann  model.  Mann  (1994)  developed  a  theory  for  the  distortion  of 
turbulence  by  a  constant  wind  shear.  The  theory  agrees  quite  well 
with  spectra  recorded  near  the  ground  on  a  windy  day.  Like  the  HGW 
model,  the  Mann  model  requires  numerical  integrations  to  compute 
most  statistics  needed  for  wave  propagation  studies. 

Models  1  to  3  actually  take  on  slightly  different  forms,  depending  on 
whether  the  turbulent  field  of  interest  is  scalar  or  vector.  Models  4  and  5 
apply  only  to  turbulent  velocity  fluctuations  (a  vector  field). 

For  each  turbulence  model,  several  statistical  functions  of  interest  must  be 
derived: 

1.  ID  correlations  and  spectra.  Although  these  are  actually  of  little  in¬ 
trinsic  interest  in  wave  propagation  calculations,  it  is  typical  of  field 
measurements  that  the  only  turbulence  data  available  are  time  series 
recorded  by  fixed  sensors  on  a  vertical  tower.  The  time  series  can  be 
converted  to  ID  spatial  series  via  Taylor's  hypothesis  (Ax  =  HAt, 
where  V  is  the  wind  speed,  and  Ax  and  At  are  the  increments  in 
space  and  time,  respectively).  From  calculations  of  spectra  and/or 
correlations  from  the  ID  spatial  series,  the  parameters  for  the  turbu¬ 
lence  model  can  be  determined.  Hence  ID  spectra  and  correlations 
are  indirectly  valuable. 

2.  2D  correlation  function.  This  function  is  needed  for  computing  the 
mutual  coherence  function,  which  describes  the  coherence  between 
the  signals  at  spatially  separated  microphones.  The  coherence  be¬ 
tween  the  microphones  is  needed  for  calculating  the  effects  of  tur¬ 
bulence  on  sensor  arrays  (Wilson,  1997a). 

3.  2D  spectrum.  Numerical  calculations  of  acoustic  propagation 
through  turbulence  are  frequently  performed  in  a  2D  plane  in  order 
to  alleviate  the  demand  on  computational  resources.  The  plane  used 
in  the  calculation  is  normally  a  vertical  one,  with  the  nominal  direc¬ 
tion  of  propagation  being  the  horizontal  coordinate.  A  2D  turbulence 
spectrum  in  the  plane  of  calculation  is  needed  to  "synthesize"  a  ran¬ 
dom  field  having  the  same  second-order  statistics  as  the  turbulence 
model. 

4.  3D  spectrum.  Most  treatments  of  volume  scattering  require  the  3D 
spectrum  for  calculations  of  the  scattered  field.  For  homogeneous 
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turbulence,  the  3D  spectrum  can  be  used  to  synthesize  a  3D  random 
field. 


5.  2D  cross  spectrum.  This  is  the  generalization  of  the  3D  spectrum 
to  vertically  inhomogeneous  turbulence.  The  2D  cross  spectrum  is 
a  spectrum  with  regard  to  the  horizontal  variables,  and  a  correla¬ 
tion  function  in  the  vertical  direction.  Inhomogeneous,  3D  random 
fields  can  be  synthesized  from  the  2D  cross  spectrum  by  a  method 
discussed  in  section  3. 

This  report  is  structured  as  follows.  Following  a  brief  exposition  of  some  of 
the  fundamental  concepts  and  equations  involved  in  modeling  turbulence 
spectra  (sect.  2),  I  discuss  a  method  for  synthesizing  a  random  turbulent 
field  from  a  prescribed  homogeneous  or  inhomogeneous  spectrum  (sect. 
3).  The  next  five  sections  each  discuss  one  of  the  five  turbulence  models 
mentioned  above.  Section  9  gives  some  examples  of  synthesized  random 
fields,  generated  by  the  various  turbulence  models. 
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2.  General  Principles 


This  section  provides  some  general  background  on  modeling  turbulence 
correlations  and  spectra. 

2.1  Correlation  and  3D  Spectral  Functions 

Let  us  define  the  correlation  function  of  a  random,  scalar  field  s  as 

R{xi,x[,X2,X2,Xs,x'3)  =  {s{xi,X2,X3)  S  {x[,X2,X3)}  ,  (1) 

where  the  angle  brackets  indicate  ensemble  averaging.  If  the  field  is  homo¬ 
geneous,  by  definition  the  correlation  function  depends  only  on  the  separa¬ 
tions  between  the  measurement  points,  Then  we  can  make  the 

simplification 

R{ri,r2,r3)  =  R  (xi,x[,X2,X2,X3,x's)  .  (2) 

If  the  field  is  also  isotropic  (i.e.,  its  statistics  are  independent  of  coordinate 
rotations),  then  the  correlation  function  of  a  scalar  depends  only  on  the 
radial  distance  r  between  the  measurement  points,  where  =  rf  r 2  -\- 
r^.  For  a  homogeneous,  isotropic  scalar  field,  one  can  therefore  define  a 
normalized  correlation  function  h  (r)  such  that 

=  i?(ri,r2,r3)  .  (3) 

The  3D  spectral  density  function  ^  (ki,  K2,  K3)  (spectrum,  for  short)  can  be 
defined  as  the  3D  Fourier  transform  of  the  correlation  function: 

^  rcsD  rcsD  pod 

<^{ki,K2,K3)  =  —  /  /  /  ii(ri,r2,r3)  (4) 

o'TT  J  —  00  J  —  00  J  —  00 

X  exp  [-i  (kiTi  -I-  K2r2  K3r3)]  dri  dr2  drs, 

/oo  noo  noo 

/  /  4>  (ki,K2,K3)  (5) 

-00  J —00  J —00 

X  exp  [i  {niri  -[-  1^2X2  «:3’’3)]  dni  dK2  dK3, 

where  Ki  is  the  wavenumber  in  the  direction  of  r*. 

Vector  quantities,  such  as  turbulent  velocity  fluctuations,  are  handled  sim¬ 
ilarly  to  scalars.  However,  the  situation  becomes  somewhat  more  compli¬ 
cated  for  vectors,  since  three  directions  are  involved  in  a  given  spectrum 
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or  correlation:  the  direction  of  the  displacement  and  the  orientations  of 
the  two  velocity  components.  (Only  the  direction  of  the  displacement  is 
involved  in  the  scalar  case.)  Let  us  define  the  correlation  function  of  a  ho¬ 
mogeneous  vector  field  as 

(n,r2,r3)  =  {ui{xi,X2,X3)uj  {x[,X2,X3))  ,  (6) 

where  the  subscripts  i  and  j  indicate  the  velocity  component.  Note  that  Rij 
is  a  tensor  quantity  having  nine  components.  We  can  define  the  spectral 
density  by  taking  Fourier  transforms  in  a  manner  analogous  to  the  scalar 
case. 

Unlike  scalar  fields,  vector  fields  (even  if  isotropic)  do  not  have  the  prop¬ 
erty  that  the  correlation  functions  depend  only  on  r.  However,  it  is  true  for 
an  isotropic  vector  field  that  Ru  {r,  0, 0)  =  R22  (0,  r,  0)  =  R33  (0,  0,  r).  Hence 
for  vectors  we  define  the  correlation  function  /  (r)  by  the  relationship 

0-2/  (r)  =  Rn  (r,  0,  0) ,  etc.  (7) 


In  this  case,  /  (r)  is  referred  to  as  the  normalized  longitudinal  correlation 
function,  so  that  it  can  be  distinguished  from  correlations  such  as  R22  {r,  0,  0) 
and  R33  (0,  r,  0),  where  the  velocity  components  and  displacement  are  per¬ 
pendicular  to  each  other.  Let  us  define  a  normalized  lateral  correlation  func¬ 
tion  g  (r)  such  that 

a^g  (r)  =  R22  {r,  0,  0) ,  etc.  (8) 

A  further  complication  when  dealing  with  turbulent  velocities,  in  compar¬ 
ison  to  scalar  fields,  is  that  the  elements  of  the  correlation  or  spectral  tensor 
cannot  be  specified  independently.  Conservation  of  mass  and  nondiver¬ 
gence  of  the  flow  impose  interrelationships  between  the  tensor  elements. 
In  fact,  the  longitudinal  and  lateral  correlations  are  related  by  the  equation 
(Batchelor,  1953) 

g  (r)  =  f  (r)  +  (9) 

The  full  correlation  tensor,  in  terms  of  the  longitudinal  and  lateral  correla¬ 
tions,  is  (Batchelor,  1953) 


Rij  {ri,r2,r3)  = 


f{r)  + 


(10) 


2.2  Energy  Spectra 

A  particularly  convenient  way  to  deal  with  some  of  the  differences  between 
the  cases  of  scalar  and  vector  random  fields  is  to  define  turbulence  models 
based  on  their  energy  spectra,  rather  than  their  correlation  functions.  The  en¬ 
ergy  spectrum  E  (k)  is  defined  such  that  its  integral  over  the  wavenumber 
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K  (where  =  Ki  +  K2  +  k^)  equals  one-half  of  the  total  variance  of  the  field. 
For  a  scalar  quantity  s, 

rco  g.2 

/  Es{K)dK=—.  (11) 

JO  ^ 

Note  from  the  Fourier  transform  relation,  equation  (5),  with  ri  =  r2  =  ra  = 

0, 

/oo  poo  poo 

/  /  ^  {ki,  K2,  K3)  dni  dK2  dK3. 

-OO  J—oo  J  —00 


Hence  we  can  find  the  isotropic  3D  spectrum  <1>  (k)  by  multiplying  Eg  {k) 
by  two,  and  then  dividing  by  the  "area"  of  a  spherical  shell  in  wavenumber 
space, 


(k) 


Eg  (k) 


(12) 


Substituting  into  equation  (5),  we  see  that  we  can  compute  the  correlation 
function  for  a  scalar  from  the  energy  spectrum  using 


h’(ri,r2,r3)  = 


/■“  Eg{K) 

Loo  27rK2 

X  exp  [i  {Kiri  +  ^2^2  -|-  dni  dK2  dn^. 


(13) 


f  —00  J  —  OO 


Because  of  isotropy,  we  can  also  set  ri  =  r,  r2  =  0,  and  ra  =  0,  before  the 
integration.  Hence 

=  /  /  TT — /  exp  (mir)  (iK2  (14) 

J-00  J-00  J-00  27tk^ 

In  the  case  of  vectors,  such  as  the  turbulent  velocities,  the  energy  spectrum 
Ey  (k)  is  defined  such  that 

fOC  ‘3^2 

/  E^{K)dK  =  —.  (15) 

Jo  2 

The  reason  for  incorporating  the  factor  3  is  that  the  total  variance  in  the  field 
equals  3a^,  if  ci^  is  defined  as  the  variance  in  just  one  of  the  three  velocity 
components.  For  a  nondivergent,  mass-conserving  flow,  the  3D  spectra  are 
related  to  the  energy  spectrum  according  to  (Batchelor,  1953) 

('^)  =  -  NKj)  •  (16) 

From  this  equation  it  can  be  shown,  analogously  to  the  scalar  case,  that 

$11  (k)  +  4>22  (A^)  +  ^33  (a^)  = 
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The  equation  for  the  correlation  function  of  a  vector  quantity,  in  terms  of 
the  energy  spectrum,  is  slightly  more  complicated  than  the  scalar  result: 


X  exp  [i  {kivi  +  K2r2  +  Ksrs)]  dni  dK2  dn^,- 

From  equation  (17),  we  have  for  vectors 

1  /■“  /■“  /■“ 


(17) 


^  /  /  Attk*  ~  '^'^2  dns,  (18) 


and 


^  ^  J  J  J  ~  ^2)  exp(iKir)  dni  dK2  d^s-  (19) 


2.3  2D  Correlation  Function 

The  2D  correlation  function  is  an  important  result  for  acoustic  and  electro¬ 
magnetic  wave  propagation  studies  (Wilson,  1998).  The  2D  correlation  func¬ 
tion  is  defined  as  the  2D  inverse  Fourier  transform  of  (0,  K2,  ^^3)/ 

/OO  POO 

/  ^  {0,K2,K3)exp[i{K2r2  + Ksrs)]  dK2dK3,  (20) 

-OO  J  —  OO 

where  -l-rl.  Equivalently,  by  Fourier  transforming  equation  (5)  with 

respect  to  ri  and  setting  ki  =  0,  one  has 

1  /■°° 

b{p)  =  —  R{ri,r2,r3)  dn.  (21) 

J  —  OO 

I  should  point  out  that  most  previous  authors  have  called  the  2D  corre¬ 
lation  function  the  transverse  correlation  function.  Unfortunately,  I  have  al¬ 
ready  used  that  term  for  the  function  g  (r).  (Taxonomically  speaking,  the 
function  g  (r)  is  a  3D  correlation  function,  with  the  spatial  separation  trans¬ 
verse  to  the  velocity  components.)  My  reason  for  adopting  the  different 
terminology  here  (besides  the  name  conflict)  is  that,  since  the  2D  spectrum 
(discussed  in  sect.  2.4)  is  the  2D  Fourier  transform  of  ii(0,r2,r3),  it  only 
makes  sense  that  the  2D  inverse  Fourier  transform  of  (0,  «;2,  ^3)  should 
be  called  the  2D  correlation  function. 

A  second  important  point  regarding  equation  (20)  is  that  the  integration 
could  just  as  well  have  been  over  the  plane  (ki,  K2)  (or  (ki,  K3),  for  that 
matter)  instead  of  (k2,  K3).  Then  the  appropriate  spectrum  in  the  integral 
would  be  4>  (ki,  K2,0),  and  the  spatial  separation  would  be  p^  =  rf  +  r^. 
The  two  definitions  are  equivalent  for  isotropic  fields. 
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The  2D  correlation  function  can  be  found  in  terms  of  the  energy  spectrum 
by  substitution  of  equation  (12)  into  (20).  The  result  is  somewhat  similar  to 
equation  (14)  for  the  normalized  correlation  function,  except  that  the  inte¬ 
gration  over  is  omitted,  and  ri  is  set  to  zero: 

b{p)  =  ^f  f  ^^-^^^exp[i{K2r2  + Ksr^)]  dK2dK3,  (22) 

J  —  OO  J—OO 

where  +  By  rewriting  this  equation  in  polar  coordinates  over  the 

plane  (^2,^3),  and  then  using  the  integral  definition  of  the  Bessel  function 
(eq  (9.1.21)  in  Abramowitz  and  Stegun,  1965),  we  can  reduce  the  double 
integral  to  a  single  one: 

b  i^p)  =  r  Jo  (^^p)  dKh.  (23) 

Jo  Kh 

The  2D  correlation  function  for  fluctuations  in  the  velocity  components  per¬ 
pendicular  to  the  displacement  can  be  found  by  substitution  of  (k),  as 
given  by  equation  (16),  into  the  definition  of  the  2D  correlation  function, 
equation  (20): 

^11  (p)  =  7“/  /  (^^2^2  +  i^srs)]  dn2  dn^  (24) 

47r  J  —  OO  J—OO 

=  I  r  dKh.  (25) 

2  JO  n-h 

Note  the  similarity  between  equations  (25)  and  (23)  —  the  only  differences 
are  the  factor  of  1/2  appearing  in  (25),  and  Ey  {nh)  replacing  Eg  (Kh)- 

The  2D  correlation  function  for  velocity  fluctuations  in  the  plane  of  the  dis¬ 
placement  is  found  by  substitution  of  either  4>22  (k)  or  4>33  (k)  into  equation 
(20): 

b±  (p)  =  ^  [  f  ^  exp  [i  {K2r2  +  ks^)]  dK2  dn^.  (26) 

47r  J  —  OO  J—OO 

In  the  above,  the  roles  of  K2  and  K3  can  be  switched  without  a  change  in 
the  result.  In  any  case,  the  function  b±  (p)  appears  to  be  of  no  importance 
in  acoustic  wave  propagation  calculations.  The  reason  is  that  it  is  primarily 
the  velocity  components  parallel  to  acoustic  wavefronts  that  affect  propa¬ 
gation.  Hence  b±  (p)  is  discussed  no  further  in  this  report. 

2.4  ID  and  2D  Spectral  Functions 

As  mentioned  in  the  introduction,  besides  the  3D  spectra,  sometimes  ID 
and  2D  spectra  are  needed  for  wave  propagation  studies.  The  ID  spectrum 
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is  simply  a  ID  Fourier  transform  of  the  correlation  function,  with  the  dis¬ 
placements  in  the  two  untransformed  directions  set  to  zero.  For  a  scalar. 


0(k) 


1 

—  /  R{r,0,0)exjp{—iKr)  dr 
ZTT  J-oc 

(J^  yoo 

—  /  /i  (r)  exp  (— iKr)  dr. 

27r  J—oc 


(27) 


I  indicate  the  normalized  ID  spectrum  for  a  scalar,  0  (k)  ,  with  the  no¬ 

tation  /r  (k). 

Of  course,  in  the  case  of  vectors,  the  relative  orientations  between  the  ve¬ 
locity  components  and  the  displacements  must  be  considered.  In  general, 
we  have 

0ii  («^i)  =  /  Rij  (ri,  0,  0)  exp  {-inin)  dri,  (28) 

^■TT  J —  oo 

and  likewise  for  the  ri  and  r2  directions.  There  are  two  cases  of  par¬ 
ticular  interest:  one  where  the  transformation  is  parallel  to  the  velocity 
components, 

f  =  TT  I  /(r)exp(-iKr)  dr  =  011  (ki) /cj^,  etc,  (29) 

J  —  OO 


and  the  other  where  it  is  perpendicular. 


Incidentally,  equations  (9),  (10),  and  (16)  can  be  used  to  show  that  the  lon¬ 
gitudinal  spectrum  is  related  to  the  energy  spectrum  according  to 


Ev  {k) 


A 

dn 


1  df{K) 
K  dn 


(31) 


A  generalization  of  the  ID  spectrum  that  is  sometimes  useful  in  wave  prop¬ 
agation  studies  is  the  ID  cross  spectrum.  Its  definition  is  the  same  as  that 
of  the  usual  ID  spectrum,  except  that  the  arguments  r2  and  rs  in  the  corre¬ 
lation  function  are  nonzero: 

0p' (Ki;r2,r3)  =  —  /  i?jj(ri,r2,r3)exp(-fKiri)  dri.  (32) 

J  — OO 

The  2D  spectrum  is  defined  as  the  2D  Fourier  transform  of  the  correlation 
function.  For  a  scalar. 


K2) 


1 

47r2 


R  {ri,r2, 0)  exp  [—i  (^iri  -|-  ^2^2)]  dri  dr2. 


(33) 
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</>ii  (ki,  H2) 


</>33(ki,K2) 


Assuming  horizontal  isotropy,  and  using  the  integral  form  of  the  Bessel 
function  (e.g.,  eq  (9.1.21)  in  Abramowitz  and  Stegun,  1965),  one  finds 

(7^  ^00 

4>{Kh)  =  7^  h{r)  Jo{Khr)rdr,  (34) 

ZTT  Jo 

in  which  kI  =  kI  +  kI  Alternatively,  from  the  definitions  of  the  Fourier 
transform  and  energy  spectrum,  it  can  be  shown  that 


4'{ki,K2)  = 


'  —00 
roc 


4>  {ki,  k.2,  K.3)  dns 
Eg  (k) 


TTK^ 


cIks- 


(35) 


The  definition  of  the  2D  spectrum  for  vectors  is  of  course  the  same  as  equa¬ 
tion  (33),  except  that  (pij  and  Rij  appear  in  place  of  (f)  and  R.  However, 
equations  (34)  and  (35)  do  not  hold.  From  equations  (9),  (10),  and  (16),  it 
can  be  shown  that 


a 

47r^  J- 


fill  A  ('  +  1) 

Ev  (k)  ^^2  _  ^2^ 


27rK^ 


a 


COO  COO 


/  /  ail 

47r  J —00  J —00 
f°°  Ey  (k)  f  2  2\ 


(36) 


(37) 


In  horizontally  isotropic  turbulence,  the  equations  for  ^22  {^1,^2)  are  the 
same  as  those  for  (pn  {ki,  K2),  except  that  the  subscripts  1  and  2  are  inter¬ 
changed.  Hence, 

</>22  (/«1,K2)  =  (pll  (k2,Ki)  . 

One  of  the  turbulence  models  considered  in  this  report  (the  HGW  model, 
sect.  7)  is  vertically  inhomogeneous.  By  definition,  this  means  that  equation 
(2)  is  invalid.  The  most  simplified  form  of  the  correlation  function  is 

R{ri,r2,z,z')  =  Rixi,xlx2,x'2,z,z)  ,  (38) 

in  which  z  =  X3  has  been  used  for  the  vertical  coordinate.  Obviously,  be¬ 
cause  of  the  inhomogeneity,  we  can  no  longer  Fourier  transform  the  cor¬ 
relation  function  with  respect  to  the  vertical  coordinate.  However,  we  can 
still  Fourier  transform  with  respect  to  the  horizontal  coordinates  (xi,  X2). 
The  resulting  spectrum  is  called  the  2D  cross  spectrum: 
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PCSD  PCO 

(j){Ki,K2]z,z')  =  j^J  J  R{ri,r2,z,z')exp[-i{Kiri  + K2r2)]  dn  dr2.  (39) 

The  2D  cross  spectrum  for  vectors  is  defined  in  the  same  way,  except  with 
subscript  ij  on  the  spectrum  and  correlation  function.  Note  that  when  2:  = 
2',  the  2D  cross  spectrum  reduces  to  the  usual  2D  spectrum.  Also,  when  the 
turbulence  is  vertically  homogeneous,  the  2D  cross  spectrum  depends  only 
on  the  separation  =  z'  —  z. 


2.5  Length  Scales 


It  is  important  that  we  be  able  to  consistently  quantify  length  scales  asso¬ 
ciated  with  given  spectra.  By  length  scale,  I  mean  a  value  associated  with 
the  size  of  eddies  from  some  part  of  the  spectrum.  One  of  the  most  useful 
length  scales  is  the  integral  length  scale,  which  characterizes  the  size  of  the 
most  energetic  eddies.  For  a  scalar  quantity,  it  is  defined  as 


C=^  R(r)  dr  = 
Jo 


h{r)  dr. 


(40) 


By  setting  k  =  0  in  equation  (27),  one  can  derive  the  following  simple  rela¬ 
tionship  between  the  integral  length  scale  and  the  ID  spectrum: 


C  =  Trh  (0)  . 


(41) 


The  integral  length  scale  and  2D  correlation  function  are  also  simply  re¬ 
lated.  From  the  Fourier  transform  relation  between  the  correlation  function 
and  the  spectrum,  and  from  equation  (12),  we  have 


1  E  (k) 

—  J  i?(ri,0,0)exp(-iKiri)  dri  =  y  J  exp  [K2r2  +  Ksrs]  dK2  dns 


Setting  r2  =  ra  =  0,  and  ki  =  0,  and  using  the  fact  that  i?(ri,  0,  0)  is  an  even 
function,  we  obtain 


1 

vr  Jo 


poo  2  /*CXD  p 

/  R{ri,0,0)dri  =  — 

Jo  ZTT  J-oc  J-00 


^  dK2  dno- 


If  this  result  is  compared  with  equations  (22)  and  (40),  it  is  now  evident  that 


TT 


C  =  ^h{t)). 


(42) 


In  the  vector  case,  the  direction  of  the  velocity  components  as  well  as  the 
direction  of  the  displacement  must  be  considered.  Two  significant  scales 
can  be  defined:  one  where  the  displacement  is  parallel  to  the  velocity. 


^11  =  /  f{r)  dr  =  TTf  (0) , 
Jo 


(43) 
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and  the  other  where  it  is  perpendicular. 


^±=  9{r)  dr  =  ng  (0) .  (44) 

Jo 

For  homogeneous,  isotropic  turbulence,  it  can  be  proven  that  (Batchelor, 
1953) 

C\l=2C±.  (45) 

Furthermore,  using  a  derivation  similar  to  the  scalar  case,  one  finds 

£||  =  46||  (0) .  (46) 


2.6  Structure  Functions  and  Parameters 

The  structure  function  for  an  isotropic,  homogeneous  scalar  field  is  defined 


as 

^(r)  =  ([5(r)-s(0)]2).  (47) 

For  a  vector,  the  structure  function  depends  on  the  direction  of  the 
displacement: 

Dij{r)  =  ([ui  (r)  -  Uj  (0)]^)  .  (48) 

From  the  definitions,  it  is  easily  shown  that 

D  (r)  =  2a^  —  2R  (r)  =  2cr^  [\  —  h  (r)]  ,  and  (49) 

Dij{r)  =  2o^  -  2Rij  (r) .  (50) 

For  vectors,  it  is  convenient  to  define 

i9||  (r)  =  2(T^  [1  —  /  (^)]  ;  arid  (51) 

D^{r)  =  2a^[l-g{r)].  (52) 


Of  particular  interest  is  the  structure  function  for  small  separations  (r  <C  £). 
Based  on  Kolmogorov's  (1941)  scaling  arguments,  the  structure  functions 
should  be  proportional  to  in  this  limit.  The  ratio  D  (r)  ,  for  small 

r,  is  called  the  scalar  structure-function  parameter  Cg.  For  vectors,  a  longitu¬ 
dinal  displacement  is  used  to  define  the  structure-function  parameter  Cf. 
Hence, 

2(7^ 

Cv  =  /  W]  •  (54) 
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Similar  to  the  ordinary  structure  function  is  the  2D  structure  function,  de¬ 
fined  for  scalars  as 


d{p)  =  2  [6(0)  -b{p)]  , 


(55) 


and  for  vectors  as 


^^11  (p) 


h  (0) 


h  (p) 


(56) 


The  2D  structure  function  is  needed  for  calculations  of  the  mutual  coher¬ 
ence  function  (see  the  introduction). 


2.7  Energy  Spectral  Approach  to  Turbulence  Modeling 

In  sections  2.1  to  2.6, 1  considered  the  interrelationships  among  various  cor¬ 
relations,  spectra,  structure  functions,  and  length  scales  that  can  be  defined 
for  random  fields.  Some  care  must  be  taken  in  devising  a  new  statistical 
turbulence  model  so  that  none  of  the  required  relationships  is  violated.  Ac¬ 
tually,  whether  one  is  dealing  with  scalar  or  vector  fields,  only  one  correla¬ 
tion  or  spectral  function  can  be  uniquely  specified;  the  others  follow  from 
the  given  relationships.  The  most  obvious  choices  for  the  specified  function 
are  the  following: 

1.  The  energy  spectrum  (k). 

2.  The  scalar  (or  longitudinal)  correlation  h  (r)  (or  /  (r)). 

3.  The  scalar  (or  longitudinal)  ID  spectral  density  h  (k)  (or  /(«))• 

The  relationship  between  E  (k)  and  /  (k)  was  given  by  equation  (31).  Of 
course,  /  (r)  and  /  (k)  are  simply  Fourier  transform  pairs.  Equations  for  all 
the  other  correlations,  spectra,  etc,  of  interest,  in  terms  of  E  (k)  and/ or  /  (r) 
have  already  been  discussed. 

In  the  following  sections  of  this  report,  I  discuss  several  turbulence  mod¬ 
els  derived  from  prescribed  forms  for  the  energy  spectrum.  The  decision 
to  use  E  (k),  instead  of  either  /  (r)  or  /  (k),  is  somewhat  arbitrary.  How¬ 
ever,  it  does  lead  to  somewhat  simpler  results  for  quantities  such  as  the  2D 
correlation  function,  an  advantage  that  I  demonstrate  shortly. 

Suppose  now  that  we  use  the  same  functional  form  for  both  the  scalar  and 
vector  energy  spectra,  by  setting* 

E^{k)  =  3Es  (k)  .  (57) 

*Of  course,  the  scalar  and  vector  fields  can  have  different  dimensions,  in  which  case 
equation  (57)  is  not  satisfied  dimensionally.  But  the  equivalence  I  imply  here  is  only  &  func¬ 
tional  equivalence  of  the  energy  spectra.  One  can  easily  overcome  any  formal  difficulties  by 
nondimensionalizing  the  fields. 
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Let  us  consider  what  additional  interrelationships  between  various  scalar 
and  vector  results  arise  from  making  this  assumption. 

With  regard  to  the  correlation  functions,  by  comparing  equations  (14),  (18), 
and  (19),  we  find 

h{r)  =  ^[f{r)+2g{r)].  (58) 

The  main  utility  of  this  equation  is  that  it  can  reduce  the  work  involved  in 
determining  the  correlation  functions.  For  example,  suppose  we  perform 
the  integrations  specified  by  equation  (18)  to  find  /  (r).  Next,  equation  (9) 
can  be  used  to  easily  determine  g  (r),  and  then  equation  (58)  can  be  used 
to  easily  determine  h  (r) .  This  procedure  is  generally  easier  than  solving 
equations  (14),  (18),  and  (19)  individually. 

By  integrating  equation  (58),  we  can  easily  derive  the  following  relation¬ 
ship  for  the  integral  length  scales: 

£  =  ^(£||+2£^)  =  ^£||.  (59) 

Fourier  transformation  of  equation  (58)  of  course  implies 

h{K)  =  ^  f  {k)  +  2g{K)  .  (60) 

A  generalized  version  of  this  equation  can  be  derived  for  the  ID  cross  spec¬ 
tra.  Using  equation  (10),  one  can  show  that  Rn  (ri,  r2,  rs)  -|-i?22  (ri,  r2,  rs)  -|- 
-^33  (?"iU2U3)  =  [/  (^)  +  2^  (r)].  It  then  follows  from  equation  (58)  that 

I?(ri,r2,r3)  =  ^  [iin  (ri,r2,r3)  -6  R22  (ri,r2,r3)  -6  i?33  (ri,r2,r3)]  .  (61) 

Fourier  transformation  of  this  equation  with  respect  to  ri  yields 

0  {tii,r2,n)  =  ^  [011  (/«iU2,r3)  +  022  (Ki,r2,r3)  -6  ©33  (Ki,r2,r3)]  .  (62) 

For  the  2D  correlation  functions,  substituting  equation  (57)  into  (24)  and 
comparing  to  equation  (22)  yields  the  following  result: 

^11  (P)  =  •  (63) 

Hence,  if  we  choose  the  same  functional  form  for  the  energy  spectrum  for 
both  the  scalars  and  vectors,  the  resulting  2D  correlation  functions  will  be 
proportional  to  one  another.  This  result  can  make  wave  propagation  calcu¬ 
lations  much  simpler. 
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If  we  apply  a  Fourier  transform  to  equation  (62)  (with  rs  =  0),  it  follows  for 
the  2D  spectra  that 

</>(ki,  K2)  =  ^  {ki,  K2)  +  </>22  {ki,  1^2)  +  hs  {ki,  K2)]  ■  (64) 


Note  also  that  </>ii  (/t2, /^i)  =  </>22  (^i,  K2)-  Equation  (64)  also  applies  to 
the  2D  cross  spectra,  after  we  make  the  notational  change  </>(Ki,/t2)  ^ 
(p{Ki,K2;z,z'),etc. 

For  the  3D  spectra,  we  have  from  equations  (12)  and  (16) 


^ij  (i^) 


(65) 


Equations  (58)  to  (65)  are  helpful  because  they  can  reduce  the  amount  of 
work  required  to  derive  various  functions  for  a  turbulence  model.  But  it 
must  be  kept  in  mind  that  they  are  valid  only  when  the  energy  spectra  sat¬ 
isfy  the  functional  equivalence  relationship,  equation  (57).  Other  modeling 
procedures  do  not  necessarily  yield  the  same  results. 
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3.  Synthesis  of  Turbulence  Using  the  Method  of  Random  Phases 


3.1  Motivation 


During  the  past  decade,  acoustic  propagation  modeling  has  now  advanced 
to  the  point  where  full-wave  solutions  can  be  computed  for  propagation 
through  a  medium  having  two-  or  three-dimensional  spatial  variability  in 
the  index-of-refraction  field  (e.g.,  Gilbert  et  al.,  1990,  Chevret  et  al.,  1996). 
What  is  still  needed,  though,  are  realistic  turbulence  fields  for  input  to  these 
models. 

Since  a  turbulence  field  cannot  be  directly  measured  to  the  fidelity  needed 
for  wave  propagation  calculations,  some  method  for  synthesizing  turbu¬ 
lence  is  needed.  Computer  simulations  of  turbulence  are  available;  in  par¬ 
ticular,  the  type  of  simulation  called  a  large-eddy  simulation  (LES)  has  been 
successfully  used  to  study  atmospheric  turbulence  for  the  past  two  decades 
(Deardorff,  1970b,  Moeng,  1984).  LES  is  regarded  as  faithful  to  the  atmo¬ 
sphere  in  most  aspects.  The  difficulty  with  LES-generated  data  is  the  reso¬ 
lution,  currently  limited  to  about  30  m  in  the  horizontal  and  10  m  in  the  ver¬ 
tical.  Supposing  that  we  need  to  resolve  the  turbulence  field  to  better  than 
1/5  acoustic  wavelength  for  reasonable  calculations,  the  LES  data  would 
be  useful  for  acoustic  propagation  modeling  only  at  acoustic  frequencies  of 
about  6  Hz  and  lower.  Therefore  alternatives  to  LES  are  needed. 

Gilbert  et  al.  (1990)  synthesized  a  planar  turbulence  field  solely  from  its  2D 
spectrum,  by  multiplying  the  square  root  of  the  spectrum  at  each  wavenum¬ 
ber  by  a  random  phase,  and  then  inverse  Eourier  transforming.  This  method 
has  also  been  used  in  other  areas  of  research.  Henceforth  I  call  this  the 
random-phase  method  (RPM).  It  is  computationally  inexpensive  and  produces 
synthetic  turbulence  having  second-order  statistics  consistent  with  the  as¬ 
sumed  turbulence  spectrum.  However,  the  fields  created  by  this  method  do 
not  necessarily  have  third-  and  higher  order  statistics  that  are  realistic  for 
turbulence.  In  particular,  RPM  cannot  capture  intermittent  features  known 
to  be  characteristic  of  turbulence.  But  if  one  is  interested  only  in  calculating 
means  and  second-order  statistics  of  the  acoustic  field,  a  method  for  pro¬ 
ducing  synthetic  turbulence  that  is  realistic  only  to  second  order  may  be 
satisfactory.  Hence  RPM  can  be  useful  in  certain  applications. 

One  shortcoming  of  RPM  as  implemented  by  Gilbert  et  al.  is  that  the  tur¬ 
bulence  model  they  used  was  an  isotropic,  homogeneous  Gaussian  model. 
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Atmospheric  turbulence,  particularly  near  the  ground,  is  known  to  violate 
all  these  properties.  Actually,  Gilbert  et  al.'s  method  can  be  used  to  synthe¬ 
size  anisotropic  non-Gaussian  turbulence.  One  would  simply  have  to  start 
with  the  appropriate  turbulence  model  for  the  spectral  density  function.  It 
is  the  inability  of  the  standard  RPM  to  incorporate  inhomogeneity  that,  in 
the  end,  limits  its  utility.  Therefore  we  need  a  method  for  generalizing  RPM 
so  that  inhomogeneous  turbulence  structure  can  be  synthesized. 

In  section  3.2,  I  describe  a  generalized  random-phase  method  (GRPM)  to 
synthesize  turbulence  for  a  given  spectral  model,  be  it  homogeneous  or 
inhomogeneous.  Next,  in  section  3.3,  I  show  how  the  method  reduces  to 
the  RPM  used  by  Gilbert  et  al.  (1990)  for  homogeneous  turbulence.  Last,  in 
section  3.4,  the  method  is  applied  to  vertically  inhomogeneous  turbulence, 
which  is  a  case  of  particular  interest  for  the  atmosphere. 

3.2  Generalized  Random-Phase  Method 

The  basic  idea  underlying  GRPM  is  to  multiply  the  empirical  orthogonal  func¬ 
tions  (EOFs)  by  random  phases  and  then  sum  them  to  form  the  synthetic 
turbulence  field.  The  reason  why  this  reduces  to  the  standard  RPM  in  the 
case  of  homogeneous  turbulence  is  that  the  EOFs  of  a  homogeneous  cor¬ 
relation  function  are  harmonic  functions  (Eumley,  1971,  Wilson,  1996).  The 
EOFs  are  defined  as  the  eigenfunctions  of  the  correlation  function: 

J  R  (x,  x')  tp  (x')  dx'  =  Xip  (x)  , 

where  x  =  (xi,  X2,  xs),  the  ip's  are  the  EOFs,  and  A  is  the  eigenvalue.  (Al¬ 
though  the  above  equation  was  specifically  written  for  the  scalar  case,  it  can 
be  used  for  vector  fields  as  well.  The  generalization  of  the  eigenvalue  prob¬ 
lem  to  vectors  and  to  multiple,  correlated  fields  is  accomplished  through 
tensors  (e.g.,  Wilson,  1996).) 

Here  is  a  step-by-step  description  of  the  GRPM: 

1.  Fourier  transform  the  eigenvalue  problem  with  respect  to  the  homo¬ 
geneous  coordinates.  Because  the  eigenvalue  equation  is  a  convolu¬ 
tion  with  respect  to  the  homogeneous  coordinates,  the  eigenvalue 
problem  is  reduced  to  a  simple  multiplication  after  the  transform,  and 
we  are  left  only  with  an  integration  over  the  inhomogeneous  coordi¬ 
nates  (Wilson,  1996).  The  result  is 

I  R  (k;  ^  (k;  ^')  d^'  =  A  (k)  ^  (k;  ^)  .  (66) 

In  the  above,  functions  that  have  been  transformed  with  respect  to 
the  homogeneous  coordinates  are  indicated  by  a  "hat,"  e.g.,  ip.  The 
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inhomogeneous  coordinates  are  indicated  by  and  the  wavenum¬ 
bers  corresponding  to  the  homogeneous  coordinates  are  indicated  by 

k. 

2.  Find  the  eigenfunctions  and  eigenvalues  in  the  transformed  domain. 
Let  us  indicate  the  solutions  as  (k;  ^),  A„  (k),  where  n  =  1,2,3,... 
The  eigenfunctions  are  orthonormal,  so  that 

J  (k;  (k;  (1$  =  bran,  (67) 

where  bmn  =  1  if  m  =  n  and  zero  otherwise. 

3.  Expand  the  transformed,  synthetic  field  s  (k;  in  terms  of  the  eigen¬ 
functions  il^n  (k;  ^): 

s  (k;  ^)  =  XI  (k)  (k;  ^)  .  (68) 

n 

Note  that  the  expansion  coefficients  an  (k)  are  random  functions. 

4.  Determine  the  coefficients  in  the  expansion  such  that 

i?(k;^,0  =  (5(k;05*(k;O)-  (69) 

(The  asterisk  indicates  complex  conjugation.) 

5.  Inverse  Fourier  transform  s  (k;  ^)  to  find  s  (x). 

The  key  to  the  procedure  is  step  4.  By  selecting  the  expansion  coefficients 
in  this  maimer,  we  assure  that  the  correlation  function  of  the  synthetic  field 
matches  the  model  correlation  function.  So,  how  do  we  choose  the  expan¬ 
sion  coefficients  such  that  equation  (69)  is  satisfied?  First  note  that  by  sub¬ 
stituting  the  expansion  (68)  into  (69),  we  have 

R  (k;  ^0  =  X  X  (k)  an  (k))  i/jm  (k;  $)  (k;  ^')  . 

m  n 

Multiplying  both  sides  by  ipj  (k;  ^'),  integrating  with  respect  to  and  in¬ 
voking  orthonormality,  we  have 

[  R  (k;  (k;  ^')  d^'  =  X  {am  (k)  a*  (k))  X  (k;  • 

Comparison  with  equation  (66)  now  yields 

X  («m  (k)  aj  (k))  X  (k;  =  Aj  (k)  (k;  ^) . 
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Finally,  multiplying  by  (k;  4),  integrating,  and  invoking  orthonormality 
one  more  time,  we  find 


(a,„(k)a;(k))  =  (5,„,,A„(k).  (70) 

This  equation  shows  that  the  expansion  coefficients  are  statistically  orthogo¬ 
nal,  and  that  the  expected  value  of  their  magnitude  squared  equals  the  cor¬ 
responding  eigenvalue.  The  simplest  way  to  choose  the  coefficients  while 
satisfying  the  orthogonality  relationship  is  to  set 

an  (k)  =  An  (k)  exp  {ijn  (k))  ,  (71) 

where  (k)  is  some  random  phase.  This  method  for  choosing  the  coeffi¬ 
cients  is  the  basis  of  GRPM. 

3.3  Homogeneous  Turbulence 

Application  of  GRPM  to  homogeneous  turbulence  is  particularly  simple. 
In  this  case,  the  Fourier  transform  is  applied  to  all  coordinates,  and  the 
transformed  eigenvalue  problem  corresponding  to  equation  (66)  is  simply 

$  {k)  Ip  (k)  =  A  (k)  ip  (k)  .  (72) 

We  make  the  eigenfunctions  orthonormal  by  setting  ip  (k)  =  1.  Hence 
<1>  (k)  =  A  (k),  and  equation  (71)  becomes 

a  (k)  =  (k)  exp  (fy  (k))  .  (73) 

The  transformed  field  is  simply  s  (k)  =  a  (k). 

3.4  Vertically  Inhomogeneous  Turbulence 

The  next  most  complicated  case  occurs  when  one  of  the  coordinates  is  inho¬ 
mogeneous.  This  is  an  important  special  case,  since  over  relatively  flat  ter¬ 
rain,  atmospheric  boundary  layer  structure  is  approximately  homogeneous 
in  the  horizontal  directions,  but  inhomogeneous  in  the  vertical.  In  this  situ¬ 
ation,  the  function  R  (k;  ^')  is  the  2D  cross  spectrum,  p  (^i,  K2',  z,  z'),  dis¬ 

cussed  in  section  2.4.  The  eigenvalue  problem  is 

j  p  {ki,  K2-,  Z,  z')  p  {ki,  K2-,  z')  dz'  =  X  {ki,  K2)  ip  {ki,  1^2]  z) .  (74) 

Generally  one  must  solve  the  eigenvalue  problem  by  discretizing  the  com¬ 
putational  domain  in  the  z-direction,  and  then  numerically  solving  the  re¬ 
sulting  matrix  equation  (Moin  and  Moser,  1989,  Wilson,  1996).  The  expan¬ 
sion  coefficients  are 
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an  {hi,k.2)  =  \J\n  («! ,  «^2)  exp  (z7„  (ki ,  ^2) )  ,  (75) 

where  n  =  1,2, . . . ,  N ,  N  being  the  number  of  grid  points  in  the  z-direction. 
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4.  Gaussian  Model 


4.1  Energy  Spectrum  and  3D  Correlations 

Starting  with  this  section,  I  consider  specific  statistical  models  for  atmo¬ 
spheric  turbulence.  Unfortunately,  there  are  no  3D  turbulence  models  that 
are  known  to  work  satisfactorily  for  a  variety  of  atmospheric  conditions. 
Hence  we  are  forced  to  consider  idealized  models. 

In  this  section  I  consider  one  such  idealization,  the  Gaussian  model.  The 
main  advantage  of  the  Gaussian  model  is  its  analytical  convenience.  But 
it  only  works  well  for  qualitatively  modeling  the  largest  scale  features  of 
the  turbulence  (the  energy  subrange).  The  von  Karman  model  (described 
in  sect.  4.2)  is  somewhat  more  general,  since  it  also  describes  the  smaller 
scale  structure  (the  inertial  subrange)  realistically. 

The  energy  spectrum,  for  the  scalar  Gaussian  model,  is  defined  as 


Es  (k) 


a'^K^L^ 


(76) 


In  accordance  with  the  discussion  in  section  2.7,  the  energy  spectrum  for 
the  vector  Gaussian  model  is  three  times  the  scalar  model: 


Ev  {k) 


(77) 


By  integration  of  equation  (31),  the  function  /  {k)  can  be  found  from  (k) 
The  result  for  the  Gaussian  spectrum  is 


The  longitudinal  correlation  function  is  found  by  calculation  of  the  inverse 
Fourier  transform: 

/(r)  =  exp  .  (79) 

The  lateral  correlation  function,  found  by  equation  (9),  is 
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Finally,  from  equation  (58),  we  have  the  scalar  correlation  function’ 


h{r)  = 


(81) 


and  its  Fourier  transform 


h  (k)  = 


L 


2  r2' 


1  + 


k^L 


2  r2' 


k^L 


exp 


Sy^r  \^  '  4  J  ^  4  ^ 

By  performing  the  integration  in  equation  (43),  we  find 


£ii  =  =  0.886L. 

2 


Hence,  from  equation  (59), 


C  =  =  0.591L. 

O 


(82) 


(83) 


(84) 


The  energy  spectrum,  as  a  function  of  kC,  is  plotted  in  figure  1. 


Figure  1.  Comparison  of 
energy  spectral  functions 
for  Gaussian,  von  Karman, 
and  Kolmogorov  models. 


Normalized  wavenumber,  k.£ 


“Previous  authors  (such  as  Ostashev,  1994)  have  used  f  {r)  =  h  (r)  =  exp  .  The 

equation  used  for  h  (r)  in  this  report  is  different  because  of  the  "energy  spectral  approach" 
I  follow.  Since  neither  equation  (79)  nor  (81)  is  known  to  agree  more  satisfactorily  with  data, 
it  is  ultimately  a  matter  of  convention  which  one  is  used. 
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4.2  ID  Cross  Spectra 

The  ID  cross  spectrum  can  be  found  by  calculation  of  the  Fourier  transform 
of  the  correlation  function  with  respect  to  ri.  The  required  integrals  are 
given  by  equations  (3.896.4)  and  (3.952.4)  in  Gradshteyn  and  Ryzhik  (1994). 
For  scalars,  one  finds 


0(«;r'2,r3)  = 


a^L 

3^ 


^  L2  + 


r 


exp 


2 

L2 


where  r2  =  +  r^.  The  results  for  vectors  are  similarly  found  to  be 

^2^  /  .2\  /  ,.2r2  .2' 


0ii(K;r2,r3)  = 
022(K;r2,r3)  = 


2^ 

a^L 

4^ 


2ro 

1 - ^  + 

L2  ^ 


2  r2' 


2  r2 


K^L 


exp 


L2 


(85) 


(86) 

(87) 


The  equation  for  ©33  {k;  r2,  r^)  is  the  same  as  the  one  for  ©22  {k]  ^2,  r3),  ex¬ 
cept  that  r2  replaces  on  the  right. 


4.3  2D  Correlation  Function 


The  scalar  2D  correlation  function  can  be  found  from  the  energy  spectrum 
by  substitution  of  equation  (76)  into  (22): 


Kp) 


exp 


4 


Jo (np)  dn. 


The  solution  to  this  integral,  from  equation  (11.4.28)  in  Abramowitz  and 
Stegun  (1965),  is 


Kp) 


a^L 


M  2 


where  M()  is  called  Rummer's  function.  From  (13.4.1)  and  (13.6.12)  in 
Abramowitz  and  Stegun  (1965),  we  see  that  M{2,l,z)  =  (1 -|- z)  x 
M  (1, 1,  z)  =  (1  +  z)  e^.  Hence  our  result  for  the  2D  correlation  function 


An  equation  equivalent  to  this  one  was  given  previously  by  Mellert  et  al. 
(1994).  The  2D  correlation  function  for  the  velocity  field  fty  (p)  is  found  sim¬ 
ply  by  multiplication  by  3/2,  as  indicated  by  equation  (63). 
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The  2D  structure  function  for  the  Gaussian  model. 


d{p)  =  2[b{p)  -6(0)] 


2cj2L 

30r 


is  plotted  in  figure  2. 


(89) 


4.4  2D  and  3D  Spectra 


The  2D  spectrum  for  a  Gaussian  scalar  can  be  found  without  much  diffi¬ 
culty  through  the  integration  in  equation  (35).  The  result  is 


2t2 


4>{Kh)  = 


a^L 

12-k 


2  r2' 


1  + 


KiL 


Similarly,  we  find  for  vectors 

^2l2 

01i(k1,K2)  = 
hsi^h)  = 


2  r2 


Stt  ^ 
IGvr 


1  + 


KfL 


exp 


exp 


k-2  r2' 


exp 


(90) 


and 


(91) 

(92) 


The  3D  spectra  are  computed  easily  from  the  energy  spectrum  (eq  (76)  and 
(77)),  with  equation  (12)  used  for  scalars  and  equation  (16)  for  vectors. 


Figure  2.  Comparison  of 
2D  structure  functions 
corresponding  to 
Gaussian,  von  Karman, 
and  Kolmogorov 
turbulence  models. 
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4.5  2D  Cross  Spectra 

One  can  find  the  2D  cross  spectrum  by  calculating  the  inverse  Fourier  trans¬ 
form  of  the  3D  spectra  with  respect  to  K3.  For  scalars,  the  resulting  integrals 
are  calculated  with  equations  (3.896.4)  and  (3.952.4)  from  Gradshteyn  and 
Ryzhik  (1994).  The  result  is 


</>(«!,  K2;r3) 


2  r2 


a^L 

127r 


1  - 


2ro 


L2 


-6 


exp 


k-2  7'2 


The  corresponding  equations  for  vectors  are 


(93) 


011  (Ki,K2;r3) 
033  {Ki,K2;r3) 


Stt 


exp 


.2  r2 


and 


(94) 

(95) 


The  equation  for  022  {ki,  ^2;  ^^3)  is  the  same  as  equation  (94),  except  that  K2 
on  the  right  side  is  replaced  by  ki. 
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5.  Von  Karman  Model 


5.1  Model  Definition  and  Energy  Spectrum 


The  von  Karman  model  is  developed  from  the  following  equation  for  the 
energy  spectrum  of  a  scalar: 


4r(i/  +  5/2) 

3\/^r(jy)  (1  +  ^2£2)!^+5/2- 


(96) 


The  parameter  £  is  a  characteristic  length  scale,  T  ()  is  the  gamma  function, 
and  V  controls  the  power-law  dependence  in  the  inertial  subrange  {k£  1). 

Generally,  we  set  =  1/3  to  obtain  Kolmogorov's  (1941)  power  law 
for  the  inertial  subrange.  Figure  1  (p  24)  compares  the  von  Karman  en¬ 
ergy  spectrum  to  the  Gaussian  energy  spectrum.  The  main  difference  be¬ 
tween  the  von  Karman  and  Gaussian  models  is  that  the  Gaussian  model 


decays  much  more  rapidly  at  large  wavenumbers.  The  vector  von  Karman 
energy  spectrum  is  simply  three  times  the  scalar  one,  so  that  equation  (15) 
is  satisfied. 


Gomplete  results  for  the  vector  version  of  the  von  Karman  model  were 
given  previously  (Wilson,  1997b).  With  these  vector  results  as  a  starting 
point,  the  equations  in  section  2.7  allow  us  to  find  the  corresponding  scalar 
results  with  little  difficulty.  For  completeness,  I  provide  here  the  important 
results  for  both  scalars  and  vectors: 


h  (r) 
h  (k) 

fir) 

/(«) 

gir) 

gi>f} 


r(5/6)  i  r  5  - 

9V^r(l/3)(l  +  ^2^2)5/6  [  l  +  ^2^2_ 


riGs)  (^) 


r(5/6)  i 

\/^r(l/3)  (1  +  k2£2)5/6’ 


riGs)  (^) 


r(5/6)  e  u  5 

V^T  (1/3)  (1  +  ^2^2)5/6  [3  -  6(1  +  k2£2) 


and 


(97) 

(98) 

(99) 

(100) 

(101) 

(102) 
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where  Ki,  is  the  modified  Bessel  function  of  the  second  kind. 
The  integral  length  scales  are 


C 

^11 


2^r  (5/6) 
3r(i/3) 

r(i/3) 


=  0.498£ 
0.747£. 


and 


(103) 

(104) 


5.2  ID  Cross  Spectra 

Equations  for  the  vector  ID  cross  spectra  were  derived  previously  (Wilson, 
1997b).  The  results  are 


011  (K;r2,r3) 
022  (K;r2,r3) 


(1/3)  [i  +  kH^) 


and 


(105) 


2a^ 


m 


5/6 


V^r(i/3)  U  + 

g-^S/el?)  ^2£2'^11/6  Vs)  +  2^^2  -^1/6(0 


(106) 


where  ^  =  (r^  +  rf)  (l  +  The  equation  for  ©33  (k;  r2,  r3)  is  the 

same  as  the  preceding,  except  that  r3  replaces  r2  on  the  right.  The  scalar 
result  follows  from  equation  (62): 


0(«:;?'2,r3)  = 


2cj2£ 

^TT^fvTT^ 


\  5/6 

9  nO  ) 


e/2 


y  ^5/6  (e)  -  Yy^2/2^ii/6  (e) 


(107) 


5.3  2D  Correlation  Function 


The  2D  correlation  function  for  the  scalar  von  Karman  model  is 


K/5) 


4c72£ 

30fr  (1/3) 


(108) 


Integrals  (3.773.6)  and  (6.726.4)  in  Gradshteyn  and  Rhyzhik  (1994)  were 
used  to  derive  this  result.  As  before,  the  vector  2D  correlation  function 
6||  (/j)  is  simply  3/2  times  the  scalar  result.  The  2D  structure  function  for  the 
scalar  von  Karman  model  is  plotted  and  compared  to  the  Gaussian  model 
in  figure  2  (p  26). 
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5.4  2D  and  3D  Spectra 


The  2D  spectra  for  the  vector  von  Karman  model,  derived  elsewhere  (Wil¬ 
son,  1997b),  are 

,  ,  4cr^K?£^ 

Dvr  (1  +  ^2 ^2) 7/3 

</>12(ki,K2)  = - ^^7!n7/3’  (^^^) 

97r  (1  +  ^t2£2)7/3 

4>13{ki,K2)  =  0.  (112) 

The  scalar  2D  spectrum,  found  with  equation  (64),  is 

n-2/'2  r  ^  ' 

The  3D  spectra  are  computed  easily  from  the  energy  spectrum,  equation 
(96),  with  equation  (12)  used  for  scalars  and  (16)  for  vectors. 

5.5  2D  Cross  Spectra 

The  2D  cross  spectra  for  the  vector  von  Karman  model  are  (Wilson,  1997b) 


(/)11  (Ki,K2;r'3)  = 


<^33(«^l,«^2;r'3)  = 


2a2^2  (4/2)4/3 

TtT  (1/3)  (1  -|-  k|£2)4/3 

[  6  2(1  -6  K^f2)  7/3  iCh) 

2a^KliHCHl2f^ 

TtT  (1/3)  (1  +  ^262)7/3  7/3  lUl  , 


where  Cft  =  (t’s/^)  y  1  +  The  equation  for  022  is  the  same  as  the  one  for 
011,  except  that  K2  replaces  ki  on  the  right.  The  scalar  2D  cross  spectrum, 
found  from  equation  (64),  is 


0(Kl,K2;r3)  = 


4(j2£2  (^^/2)4/3 

37rr(l/3)(l  +  4£2)4/3 

X  y  IT4/3  (4)  -  ^  ^2^2)  1^7/3  (Ch) 
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6.  Kolmogorov  Model 


6.1  Energy  Spectrum  and  Model  Definition 

The  Kolmogorov  model  is  based  on  scaling  hypotheses  that  apply  only  to 
the  inertial  subrange:  that  is,  for  p  C  in  the  spatial  domain,  or  »  1 
in  the  wavenumber  domain.  Actually,  the  von  Karman  model  discussed 
in  the  previous  section  has  an  inertial  subrange  satisfying  Kolmogorov's 
hypotheses.  Therefore  the  approach  I  take  for  deriving  the  Kolmogorov 
model  is  simply  to  calculate  limiting  forms  of  the  von  Karman  model.  It 
needs  to  be  pointed  out,  though,  that  this  procedure  produces  a  possible 
Kolmogorov  model;  there  is  no  unique  Kolmogorov  model.  Kolmogorov's 
hypotheses  provide  only  proportionality  relations;  the  numerical  constants 
appearing  in  the  equations  in  this  section  depend  in  part  on  assumptions 
intrinsic  to  the  von  Karman  model. 

The  energy  spectrum  for  scalars  is  the  large  KAlimit  of  equation  (96): 


Eg  (k) 


4r(17/6) 

3V^r(i/3) 


a^£  (k£) 


-5/3 


(117) 


(Since  £  and  C  are  the  same  order,  »  1  is  equivalent  to  »  1.)  This 
scalar  energy  spectrum  is  plotted  in  figure  1.  As  before,  we  choose  Ey{K)  = 
3Es  (k).  In  the  next  section,  I  derive  expressions  for  the  structure-function 
parameters  (defined  later  in  eq  (130)  and  (131)).  These  enable  us  to  write 
the  energy  spectra  as 


Eg  {k) 


Ey  (k) 


~  0.2074C;V-V3,  and  (118) 


Equations  (118)  and  (119)  (but  not  eq  (117))  are  the  ones  that  would  normally 
be  used  for  the  Kolmogorov  model.  The  reason  is  that  equation  (117)  ex¬ 
plicitly  contains  ci^  and  £,  two  parameters  that  depend  primarily  on  the 
large-scale  structure  of  the  flow.  On  the  other  hand,  equations  (118)  and 
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(119)  contain  structure-function  parameters,  which  describe  the  small-scale 
(inertial  subrange)  structure  of  the  flow.* 

For  the  ID  spectra,  we  find 


In  terms  of  the  structure-function  parameter,  the  ID  scalar  spectrum  is 

(k)  =  c.  0.1244C2«-5/3.  (122) 

The  equation  for  the  ID  longitudinal  vector  spectrum  is  the  same  as  equa¬ 
tion  (122),  except  with  replacing  (7^.  The  ID  transverse  spectrum  is  4/3 
times  the  longitudinal  spectrum. 


*The  Kolmogorov  spectral  equations  given  by  Tatarskii  (1971)  and  Ostashev  (1994)  differ 
somewhaf  from  equations  (118)  and  (119).  For  example,  the  3D  spectral  density  of  a  scalar 
is  given  as 


4>(k) 


5V3r  (2/3)  2  -11/3 


(a) 


On  fhe  other  hand,  since  I?  (k)  =  I?®  (k)  1‘1-kk^ ,  equation  (118)  implies 


4>(^) 


2^/^  r(17/6)  2  -11/3 
ll7r3/2  r(2/3)  " 


(b) 


To  prove  that  this  equation  is  equivalent  to  the  one  given  by  Tatarskii  and  Ostashev,  we  first 
use  the  recursion  property  of  the  gamma  function,  T  (2  -I-  1)  =  ^T  (2).  Hence, 


r 


From  the  duplication  formula  for  fhe  gamma  function  (eq  (6.1.18)  in  Abramowitz  &  Stegun, 
1965),  we  have 


Furthermore,  from  the  reflection  formula  (eq  (6.1.17)  in  Abramowifz  &  Stegun,  1965), 


Putting  these  various  results  together,  we  find 


r(17/6)_55  V3  /2\ 
r(2/3)  36  22/3^  Vs/- 

Tatarskii  and  Ostashev's  equation  (a)  follows  immediately  upon  substitution  of  fhis  result 
into  equation  (b). 
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To  find  the  limiting  forms  of  the  correlations  for  small  v  jt,  the  following 
expansion  of  the  Bessel  functions  for  small  arguments  can  be  used: 


2 


2v  V2/ 


(123) 


Note  that  /  (r),  g  (r),  and  h  (r)  in  the  von  Karman  model  all  have  the  form 


ip{r) 


2 

r(i/3) 


(124) 


where  a  =  0  for  /  (r),  a  =  1/2  for  g  (r),  and  a  =  1/3  for  h  (r).  For  the  inertial 
subrange,  then,  the  general  result  is 


ip{r) 


1  -  (3  +  2a) 


r(2/3) 

r(i/3) 


(125) 


6.2  Structure  Functions  and  Parameters 


The  structure  functions  are  found  by  substitution  of  equation  (125)  into 
equations  (49),  (51),  and  (52).  The  results  are 


D{r) 
Du  (r) 


D±  (r) 


22r(2/3)  / 

'  r  \ 

2/3 

(126) 

3r(i/3)  [ 

6a^T{2/3) 

/  j. 

\2/3 

j  ,  and 

(127) 

r(i/3) 

V^. 

8a2r(2/3) 

r(i/3) 

/  j. 

{Ye 

^2/3 

(128) 

Using  equations  (10)  and  (50),  we  have  more  generally  for  the  velocity 
fluctuations 


(^) 


2a2r(2/3)  /  r  \2/3  / 

r(l/3)  \ 


(129) 


The  structure-function  parameter  is  defined  as  the  ratio  D{r) / r^/3  for  small 
r.  Hence  the  structure-function  parameter  for  scalars  is 


22cj2r(2/3) 

3r(l/3) 


~  2.335cJ2r2/^ 


(130) 


and  for  vectors 


Cf, 


6a2r(2/3) 

r(i/3) 


~  i.giicj^r^/^. 


(131) 
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Therefore  equation  (129)  can  be  rewritten  as 


fin 


(r)  =  (^4  —  cos^ 


=  Czr 


2^2/3 


1  +  (1/3)  sin^  a 


where  cos  a  =  r^/r. 

For  small  separations,  the  scalar  2D  correlation  function  is 


b{p) 


2cj2£r(5/6)  r  nr (1/6) 
3V^r(i/3)  “  yr(5/6)  \^) 


The  2D  structure  function  is  hence 

Ua^iT  {1/6)  ( 


2[h  {6) -b{p)]  = 


For  vectors. 


15  V^r(i/3) 


^  2.143a2r2/3p5/3. 


d|l  {p)  =  '^d{p)  ~  3.214cj2r2/3p5/3^ 


In  terms  of  the  structure-function  parameters,  we  can  write 


d{p) 
d\\  {p) 


(132) 

(133) 

(134) 

(135) 

(136) 

(137) 


Equation  (136)  is  plotted  in  figure  2  (p  26). 


7.  HGW  Model 


The  HGW  model  is  essentially  a  von  Karman  model,  modified  to  account 
for  the  blocking  effect  of  the  ground  on  the  turbulence.  As  a  result  of  ground 
blocking,  the  vertical  velocity  of  an  eddy  must  vanish  at  ground  level.  This 
causes  the  turbulent  velocity  statistics  to  be  statistically  inhomogeneous 
in  the  vertical  direction.  Since  I  describe  the  HGW  model  in  detail  else¬ 
where  (Wilson,  1997b),  here  I  only  summarize  the  main  results,  and  show 
how  they  can  be  used  for  wave  scattering  calculations  and  random  field 
synthesis. 

In  previous  sections,  orientation  of  the  coordinate  system  was  arbitrary. 
Because  the  HGW  model  is  anisotropic,  we  must  be  specific  in  this  section 
by  using  X3  (that  is,  z)  for  the  vertical  axis. 

The  basis  of  the  HGW  model  is  the  following  equation  for  the  2D  cross 
spectrum  of  the  turbulence: 

(l)ij  {ki,K2]Z,z')  =  {ki,K2-,\z'  -  z\)  (138) 

-h  (ki,  K2)  (piP*  {ki,  K2;  z) 

+  (ki,  K2)  (pip  {ki,  K2',  z') 

+  (ki,  K2)  rrij  {ki,  K2)  (jpP  («:i,  K2;  0) , 


where 


nij  (ki,  K2) 


j  =  1  or  j  =  2 

-1,  j  =  3 


(139) 


The  superscripted  "{H)”  in  the  equations  above  means  the  homogeneous 
spectrum:  that  is,  the  spectrum  that  would  be  measured  if  the  boundary 
(ground)  were  not  present.  When  the  third  argument  of  (plj  ^  (ki,  K2H3)  is 
zero,  the  2D  cross  spectrum  becomes  the  ordinary  2D  spectrum,  so  that 
equations  (109)  to  (112)  can  be  used.  For  nonzero  argument,  the  following 
results  from  Wilson  (1997b)  are  needed  in  addition  to  equations  (114)  and 
(115): 


(piP  {Ki,K2;r3)  = 
(p[P  {Ki,K2;r3)  = 


^r{i/3){i  +  4Pg^ 


2ia^Ki£^  (a/2)^/^ 


7rr(l/3)  (1  +  4^2) 


11/6  ■ 


(140) 

(141) 
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The  main  difficulty  in  applying  the  HGW  model  is  that  equation  (138)  can¬ 
not  be  Fourier  transformed  analytically  One  is  therefore  forced  to  use  nu¬ 
merical  methods. 

7.1  ID  Cross  Spectra 

The  ID  cross  spectra  must  be  calculated  numerically  They  are  very  useful 
for  synthesizing  turbulence  in  a  vertical  plane.  As  an  example,  suppose 
we  wish  to  synthesize  the  ui  velocity  component  in  the  (xi,  X3)-plane.  The 
EOFs  required  by  the  GRPM  are  the  eigenfunctions  of  the  3D  correlation 
function  Ru  {ri,0,  z,  z').  We  can  calculate  this  function  by  taking  the  2D 
inverse  Fourier  transform  of  equation  (138)  with  i  =  j  =  1,  and  evaluate 
it  at  r2  =  0.  The  next  step  is  to  Fourier  transform  the  correlation  function 
with  respect  to  the  homogeneous  coordinates,  in  this  case  ri .  The  resulting 
function  is  the  ID  cross  spectrum  ©n  {ki;  0,  z,  z').  It  is  easily  shown  that 

<dn{Ki]0,z,z')  =  —  Rn{ri,0,z,z')exp{-iKiri)  dri 
■^TT  J  —  OO 

/CO 

(j)ii{Ki,K2;z,z')  dK2.  (142) 

-OO 


The  second  form  is  more  useful,  since  it  allows  ©n  to  be  determined 
from  a  simple  ID  integration  of  (^n,  which  is  known  from  equation  (138). 
Hence  for  synthesis  in  a  vertical  plane,  the  ID  cross  spectral  function 
©11  (ki;  0,  z,  z')  plays  the  role  of  R  in  the  procedure  described  in  section  3.2. 
Furthermore,  k  ^  ki,  and  $  ^  z.  Once  these  identifications  are  made,  it  is 
straightforward  to  apply  the  GRPM. 

7.2  2D  Correlation  Function 

For  homogeneous  turbulence,  the  2D  correlation  function  can  be  deter¬ 
mined  from  either  equation  (20)  or  (21).  Because  the  3D  spectrum  in  equa¬ 
tion  (20)  is  undefined  in  inhomogeneous  turbulence,  however,  equation 
(21)  is  the  only  definition  directly  compatible  with  the  HGW  model. 

Let  us  consider  the  calculation  of  fty,  for  which  the  integration  and  velocity 
components  are  aligned.  Using  the  xi-axis  for  the  direction  of  propagation, 
we  have 

1 

b\\  {r2,z,z)  =  —  Rn  {ri,r2,z,z')  dri.  (143) 

J— OO 

Equivalently,  from  the  Fourier  transform  definition,  equation  (4),  one  has 

poo 

b\\{r2,z,z')  =  /  (j)ii{0,K2-,z,z')exp{iK2r2)dK2.  (144) 
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Because  mi  (0,  K2)  =  0,  only  the  first  term  in  equation  (138),  represent¬ 
ing  the  contribution  to  the  spectrum  that  is  unaffected  by  the  ground,  is 
nonzero.  Therefore  the  2D  correlation  function  for  the  HGW  model  is  the 
same  as  the  von  Karman  model,  equation  (108). 

7.3  2D  and  3D  Spectra 

It  is  a  simple  matter  to  obtain  the  2D  spectrum  in  a  horizontal  plane  for  the 
HGW  model:  one  just  sets  z  =  z'  in  equation  (138).  On  the  other  hand,  the 
2D  spectra  in  vertical  planes  and  the  3D  spectrum  are  undefined  because 
of  vertical  inhomogeneity. 
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8.  Mann  Model 


The  Mann  model  is  the  result  of  applying  rapid  distortion  theory  to  a  con¬ 
stant  shear  layer.  The  main  idea  underlying  rapid  distortion  theory  is  that 
the  shear  distortions  to  an  eddy  occurring  over  short  time  intervals  can  be 
modeled  by  linearized  Navier-Stokes  equations.  Although  a  full  discussion 
of  rapid  distortion  theory  and  Mann's  model  is  beyond  the  scope  of  this  re¬ 
port,  the  relevant  equations  are  summarized  here. 

Mann's  equations  for  the  3D  autospectra  in  a  uniform  (constant  gradient) 
shear  layer  are 

^11  (k)  =  [kl  -kj-  2k,ksoCi  +  {kl  +  fci)  Cf]  ,  (145) 

‘I>22  (k)  =  [kl  -  kl  -  2k2k3oC2  +  {kl  +  kl)  Cl]  ,  and  (146) 

In  the  equations  above,  E  (ko)  is  the  initial  (before  the  onset  of  shear  dis¬ 
tortion)  energy  spectrum.  The  isotropic  vector  von  Karman  energy  spec¬ 
trum  is  used  for  E  (ko).  The  initial  wavenumber  is  ko  =  (ki,  k2,  kzo),  where 
^30  =  ks  —  f3ki,  and  (3  is  called  the  nondimensional  eddy  lifetime.  It  is  given  by 
(Mann,  1994,  Wilson,  1998) 


The  (i  are  given  by  the  equations 

Cl  =  Ci-^C2,  C2  =  ^Ci  +  C2, 

kl  kl 

where 


_  (3kl  {kl  -  2klQ  +  (3kik-io) 

k^  {kl  +  kl) 

k2kl  \ Pki  {kl  +  kl)^^'^ 

{kl  +  kl)^^  [  kl-pk3oki 


(149) 

(150) 

(151) 
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9.  Examples  of  Synthesized  Random  Turbulence 


I  provide  here  some  examples  of  synthesized  turbulence  fields,  generated 
by  the  various  turbulence  models  discussed  in  previous  sections. 

9.1  Large-Scale  (Buoyantly  Driven)  Turbulence 

The  turbulence  parameters  chosen  for  this  set  of  examples  are  intended 
to  be  characteristic  of  the  large  eddies  generated  by  buoyancy  instabili¬ 
ties,  such  as  occur  when  the  sun  heats  the  ground  during  the  day  Specif¬ 
ically,  the  parameter  values  for  the  HGW  model  are  those  suggested  by 
Wilson  (1997b):  =  O.SStc^,  where  wl  is  called  the  mixed-layer  velocity  scale 

(Deardorff,  1970a),  and  £  =  0.232;*,  where  Zi  is  the  boundary-layer  inversion 
height.  Besides  the  HGW  model,  the  models  used  to  make  the  synthesized 
fields  for  this  case  are  the  homogeneous  von  Karman  (sect.  5)  and  Gaus¬ 
sian  (sect.  4)  models.  The  parameters  used  in  the  von  Karman  model  are 
the  same  as  in  the  HGW  model.  For  the  Gaussian  model,  cr^  =  0.35rt;J  and 
L  =  2r  (5/6)  /T  (1/3)  I  =  0.1942;*  were  used.  The  reason  for  choosing  this 
value  of  L  is  that  it  results  in  matching  integral  length  scales  for  the  von 
Karman  and  Gaussian  models,  as  is  evident  from  a  comparison  of  equa¬ 
tions  (83)  and  (104). 

Figures  3  to  5  show  the  Gaussian,  von  Karman,  and  HGW  models,  respec¬ 
tively.  All  the  figures  depict  vertical  planes.  The  field  synthesized  is  the 
in-plane  horizontal  velocity  component.  The  resolution  of  the  simulated 
fields  is  150  grid  points  in  the  vertical  direction  and  200  in  the  horizontal. 
The  Gaussian  field  was  computed  in  just  a  few  minutes  on  a  Sun  Ultra- 
Sparc  workstation;  the  von  Karman  field  required  several  hours,  mainly 
because  of  the  Bessel  function  evaluations;  and  the  HGW  field  required 
several  days,  because  of  the  numerical  integration  to  determine  the  ID 
cross  spectrum  from  the  2D  cross  spectrum.  In  general,  most  of  the  com¬ 
putation  time  involves  determining  the  ID  cross  spectra;  the  eigenanalysis 
to  determine  the  EOFs  is  fast,  as  is  application  of  the  GRPM.  Because  all 
three  models  have  axial  symmetry,  the  azimuthal  orientation  of  the  vertical 
plane  relative  to  the  wind  direction  plays  no  role.  The  ID  cross  spectra  (eq 
(86)  and  (105)  with  r2  =  0,  and  eq  (142))  were  used  to  synthesize  the  fields. 
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Figure  3.  Synthesized 
turbulent  velocity  field  for 
buoyanfly  driven 
turbulence,  generated  by  a 
homogeneous  Gaussian 
model. 


*(m) 


Figure  4.  Synthesized 
turbulent  velocity  field  for 
buoyanfly  driven 
turbulence,  generated  by  a 
homogeneous  von  Karman 
model. 
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Figure  5.  Synthesized 
turbulent  velocity  field  for 
buoyanfly  driven 
turbulence,  generated  by 
Hunt/ Graham/Wilson 
model. 
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The  Gaussian  model  produces  a  much  more  smoothly  varying  field  than 
the  HGW  and  von  Karman  models.  Such  smoothness  is  unrealistic  for  high- 
Reynolds-number  turbulence  such  as  the  atmosphere.  The  HGW  and  von 
Karman  models  are  much  more  realistic,  and  yield  quite  similar  results. 
There  is  a  slight  enhancement  of  the  velocity  fluctuations  near  the  ground 
in  the  HGW  model,  which  results  from  the  ground  blocking  effect.  (The  dif¬ 
ferences  between  the  HGW  and  von  Karman  models  are  much  more  dra¬ 
matic  for  the  vertical  velocity  component.  The  horizontal  velocity  was  plot¬ 
ted  here,  since  it  is  more  important  for  horizontal  acoustic  propagation.) 


9.2  Small-Scale  (Shear-Driven)  Turbulence 

The  parameters  values  chosen  for  this  set  of  examples  are  characteristic 
of  turbulence  generated  by  surface-layer  wind  shear.  The  integral  length 
scale  for  this  type  of  turbulence  is  proportional  to  the  height.  Five  example 
synthesized  turbulence  fields  are  shown  in  figures  6  to  10.  As  before,  the 
synthesized  field  is  the  in-plane  horizontal  velocity  component.  The  reso¬ 
lution  of  the  simulated  fields  is  200  grid  points  in  the  vertical  direction  and 
300  in  the  horizontal. 

The  first  of  the  synthesized  fields  (fig.  6)  was  created  with  an  inhomoge¬ 
neous  Gaussian  model  with  =  2.97u^,  where  is  called  the  friction 
velocity,  and  L  =  1.342;,  where  z  is  the  height.  The  second  synthesized  field 
(fig.  7)  is  an  inhomogeneous  von  Karman  model,  with  =  2.97u^  and 
£  =  l.GOz.  (The  parameter  values  for  both  models  are  the  ones  derived  for  a 
shear-driven  surface  layer  by  Wilson,  1998.)  The  inhomogeneous  Gaussian 
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Figure  6.  Synthesized 
turbulent  velocity  field  for 
shear-driven  turbulence, 
generated  by  an 
inhomogeneous  Gaussian 
model  with 

height-dependent  length 
scale. 


x{m) 


Figure  7.  Synthesized 
turbulent  velocity  field  for 
shear-driven  turbulence, 
generated  by  an 
inhomogeneous  von 
Karman  model  with 
height-dependent  length 
scale. 


x(m) 


and  von  Karman  fields  both  clearly  show  small  structures  near  the  surface 
that  gradually  merge  into  larger  ones  away  from  the  surface.  However,  the 
Gaussian  model  produces  an  unrealistically  smooth  field. 
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Figure  8  is  a  second  example  of  a  Gaussian  field,  differing  from  figure  6 
in  that  the  length  scale  has  been  assigned  a  constant  (height-independent) 
value  of  L  =  1.1m.  Although  this  value  for  L  is  rather  unrealistic,  it  has 
often  been  used  previously  for  acoustic  wave  scattering  studies  (Daigle  et 
al.,  1983,  Daigle  et  al.,  1986,  Gilbert  et  al.,  1990,  Ghevret  et  al.,  1996).  The  re¬ 
sulting  homogeneous  field  lacks  the  realistic  merging  of  smaller  structures 
into  larger  ones  that  was  evident  in  figures  6  and  7. 

The  final  two  s}mthesized  fields  I  consider,  figures  9  and  10,  were  created 
with  the  Mann  model.  For  reasons  discussed  in  a  previous  report  (Wilson, 
1998),  the  parameter  values  used  in  the  Mann  model  are  T  =  3.53,  cr^  = 
1.88rt^,  and  £  =  O.SOSz.  Recall  from  section  8  that  the  spectra  in  the  Mann 
model  depend  on  the  orientation  relative  to  the  wind.  Figure  9  shows  an 
along-wind  vertical  plane,  while  figure  10  is  a  crosswind  vertical  plane. 
Qualitatively,  both  fields  are  quite  similar  to  the  von  Karman  model  (fig. 
7).  Although  it  is  difficult  to  discern  visually,  there  is  a  statistical  tendency 
for  larger  structures  in  the  along-wind  direction  of  the  Mann  model  than 
in  the  von  Karman  model.  Some  decrease  of  the  variance  in  the  crosswind 
direction  is  visually  evident.  It  is  interesting  that  the  computations  for  the 
Mann  model  required  about  one  day,  which  is  significantly  less  than  the 
HGW  model,  even  though  both  models  require  numerical  integrations.  The 
Mann  model  ends  up  being  faster  because  it  is  based  mostly  on  simple 
algebraic  functions,  whereas  the  HGW  model  contains  Bessel  functions. 


Figure  8.  Synthesized 
turbulent  velocity  field  for 
shear-driven  turbulence, 
generated  by  a 
homogeneous  Gaussian 
model  with  a  constant 
(height-independent) 
length  scale. 
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Figure  9.  Synthesized 
turbulent  velocity  field  for 
shear-driven  turbulence, 
generated  by  Mann 
rapid-distortion  model. 
Shown  is  a  vertical  plane 
parallel  to  wind; 
synthesized  field  is 
horizontal  wind 
fluctuation  parallel  to 
mean  wind. 
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6 
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Figure  10.  Synthesized 
turbulent  velocity  field  for 
shear-driven  turbulence, 
generated  by  Mann's 
rapid-distortion  model. 
Shown  is  a  vertical  plane 
perpendicular  to  wind; 
synthesized  field  is 
horizontal  wind 
fluctuation  perpendicular 
to  mean  wind. 
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10.  Concluding  Remarks 


A  comprehensive  set  of  equations  for  second-order  statistical  models  of 
turbulence  is  derived  in  this  report.  Future  research  will  involve  application 
of  these  results  to  wave  propagation  calculations. 

The  generalized  random-phase  method  (GRPM)  described  in  this  report  for 
synthesizing  inhomogeneous  random  fields  provides  much  more  realistic 
turbulence  fields  for  use  in  wave  propagation  calculations  than  do  previous 
methods.  The  method  could  also  be  applied  to  turbulent  transport  and  dif¬ 
fusion  calculations.  When  GRPM  is  used  with  the  HGW  and  Mann  turbu¬ 
lence  models,  the  resulting  synthesized  fields  have  second-order  statistics 
very  close  to  the  known  properties  of  atmospheric  turbulence. 
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Acronyms 


ID  (2D,  3D) 

CBL 

EOF 

GRPM 

HGW 

RPM 


one-  (two-,  three-)  dimensional 
convective  boundary  layer 
empirical  orthogonal  function 
generalized  random-phase  method 
Hunt/ Graham/Wilson  model 
random-phase  method 
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